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M accurac:y study of central finite difference methods 


IN SECOND ORDER BOUNDARY VALUE PROBLEMS 

By- 

Nancy Jane Cyrus 
ABSTRACT 

An accuracy study is made of central finite difference methods 
for solving boundary value problems which are governed by second 
order differential equations with variable coefficients leading to 
odd order derivatives. Three methods are studied through applications 
to selected problems. Definitive expressions for the error in each 
method are obtained by using Taylor series to derive the differential 
equations which exactly represent the finite difference approximations. 
The resulting differential equations are accurately solved by a 
perturbation technique which yields the error directly. A half 
station method, which corresponds to making finite difference 
approximations before expanding derivatives of function products in the 
differential equations, was found superior to two whole station methods 
which correspond to expanding such products first. 
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IV, INTRODUCTION 

In the mathanatl cal analysis of many physical houndarj'- value, 
problems, such as beams, plates and she31.s in structural analysis, 
the governing differential equations axe often solved by approximating 
the derivatives by finite differences and solving the resulting system 
of algebraic equations on a digital computer. In the analysis of 
complicated structures the number of simultaneous equations resulting 
from finite differences may be large enough to exceed the capacity of 
the computer or to introduce round-off error in obtaining a numerical 
solution. For such problems, it is important to keep the n\imber of 
algebraic equations at a minimum and the accuracy of the difference 
procedure can be a critical item in obtaining meaningful results. In 
reference 2, for example, it ms found that accurate answers for the 
stresses in a shell structure could not be obtained by using certain 
finite difference approximations unless the mesh spacing was smaller 
than machine capacity permitted. 

The most popular difference approximations used in boundary value 
problems are the central difference approximations which are given in 
textbooks on numerical methods. There are alternate formulations of 
central differences which can be used when odd order derivatives occur 
in the differential equation and these alternate formulations give 
different answers. It was shown in reference l4 that for a circular 
plate symmetrically loaded, approximating the differential equation 
by central differences led to a nonsymmetric matrix instead of the 
expected symmetric matrix. Furthermore, the answers in no way 


t 


. 6 - 

resembled the known solutions to the problem and the central difference 
equation was singulsir at the center of the plate, a physically real 
point in the problem. 

The purpose of this paper is to investigate the accuracy of the 
three alternate forms of central finite difference approximations as 
applied to boundary value problems. An approach for studying the 
accuracy of finite difference methods is presented and utilized. The 
study is confined to linear second order boiindary value problems of a 
certain type but the approach and conclusions are applicable to a wide 
class of boundary value problems. 
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V. General Discussion of Error 
Types of Error 

The use of finite difference approximation formulas to obtain 
numerical solutions to differential equations leads to errors which 
can be classified as three types: (l) round-off error, (2) inherited 

error, and (3) truncation or discretization error. Round-off error 
is a calculation error resulting from using a finite number specified 
by n correct digits to approximate a number which requires more than 
n digits for its exact specification. Roimd-off error increases with 
the number of calculations required to get an suiswer. The inherited 
error is the contribution to the error due to the total error at a 
preceding step. This may result from using a step-by-step procediore 
in which each step uses the restilt from the prervlous step. 

Truncation error, or discretization error as it is sometimes 
called, comes from approximating or replacing the continuous problem 
by a discrete model. Discretization error is decreased by using 
smaller increments; but as increment size decreases, the number of 
steps talcen increases, calculations increase, and the danger that 
round-off error will build up to substantial proportions grows. In 
any problem that is short enough to permit hand computation, it is 
usually possible to carry enough places so that round-off error can 
be neglected. In extended computations using computing machines 
round-off error can be serious. 

All three types of error can occur when a boundary value differ- 
ential equation is solved by reducing it to an initial value problem 
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and then solving it by one of the step-by-step procediires for initial 
value problems (ref, 10) . If a boundary value problem is solved by 
replacing the differential equation by central difference equations, 
taking into account the boundary conditions at both ends, and thus 
obtaining a set of simultaneous algebraic equations, irlierited error 
does not exist as a separate entity. For such problems round-off and 
truncation error are the only separable effects. Round-off error, 
while it can be important in a practical problem utilizing large 
numbers of simultaneous equations, is not considered here. 

Literatiire Survey 

Nijmerous studies have been reported in the literature dealing 
with errors resulting from the use of numerical methods to approximate 
the solutions to linear and nonlinear ordinary and partial differential 
equations governing boixndary value problems. A common way to solve a 
boundary value problem approximately is to reformulate the problem as 
an initial value problem and solve it using numerical integration. 
Consequently most of the error studies in the literature deal with 
initial value problems. However, some comments on a few important 
papers and books which do treat errors in boundary value problems are 
given here. 

CoUatz (ref, 3) gives methods for solving boundary value 
problems directly and for obtaining estimates of the discretization 
error. This is acconqilished by first expanding the difference 
equations in Taylor series, then deriving a system of equations for 
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the errors, estimating higher order derivatives in some way and then 
solving the system of error equations to obtain error bounds. 

In Modern Computational Methods, reference 15, a difference 
correction is added to the central difference approximations for the 
derivatives. A first approximation solution to the resulting system is 
obtained by neglecting the difference correction and solving the 
resulting algebraic equations. Then considering the difference 
correction a successive correction method is used to obtain corrections 
to the first approximation solution. The process is continued until 
there is no change in the numerical solution. 

Many methods of error analysis of boundary value problems in 
partial differential equations are also applicable to ordinary differ- 
ential equations. In the classic method developed by Gerschgorin 
(ref. 6), the discretization error is estimated by the use of a special 
method which he calls the majorant method. This method is also 
discussed by Collatz (ref. 5) and Forsythe and Wasow (ref. 5) • 

Roudebush (ref. 15) uses an error analysis of the Gerschgorin type to 
show that the order of discretization error in ordinary differential 
equations and parabolic and elliptic partial differential equations is 
unaffected by a finite number of discontinuities in the coefficients 
of the differential equation. In this paper he derives some higher 
order finite difference approximations and shows that when these 
approximations are used the order of the discretization error is 


improved . 
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Bramble and Hubbard (ref. l) have included the work of 
Gerschgorin (ref. 6) and Collatz (ref. 3) as special cases in their 
theorem for estimating error in the Dirchlet problem for elliptic 
equations . 

In many studies of physical problems approximate methods are 
judged with the knowledge of what the correct solution should be. In 
Chuang and Veletsos (ref. 2), for example, two finite difference 
methods are used to obtain approximate solutions to the partial 
differential equations governing the deformation of cylindrical shell 
structures. One method gives results which are unacceptable, even as 
design data, while the other method gives a satisfactory solution. 

Round-off error resulting from the solution of tridiagonal 
matrices, which result from the use of central difference methods 
in some boundary value problems, is not the concern in the present 
paper but has been treated to some extent in the literature. Von 
Neumann and Goldstine (ref. 19) establish an error bounds for which 
solutions by the elimination method is valid. Turing (ref. l8) 
discusses different matrix methods and gives round-off errors for the 
Jordan, Gauss and Choleski methods. Wilkinson (ref. 20) also gives 
estimates of round-off error in matrix solutions, while Lowan (ref. 11) 
deals specifically with tridiagonal matrices. 
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VI . DEVELOPMENT OF FINITE DIFFEEEKCE OPERATORS 

Finite difference operators can be obtained by several methods. 
Three common procedures are given and used to derive the difference 
approximations investigated in the study. 

Polynomial Approximation 

One method of obtaining approximate values for the derivatives of 
a fxinction which is known at a discrete number of points consists of 
fitting the given points with an appropriate polynomial, whose 
derivatives are then obtained. Referring to figure 1 the problem is 
to find the derivatives of the function which passes through the 
given points Jx^, y^y ^x^, yj ' * * Values of the 

function are known at these points or stations. 

Lagrange's interpolation formula can be specialized to fit a 
polynomial through a certain number of points . Let there be given 
values of the ordinates y , y, * • • y of the function y - f(x) 
at the (n + l) points x^, x^, • • • x^. The polynomial of the 
nth degree through these points may be written in the form 
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The equation for a polynomial passing through three points 
separated by equal increments h and with the origin at x = is 
obtained from equation ( 6 . 1 ) 


- i'c i; * S 


=^2) * 


2 h 


(y„ - 2^1 + y^) 


(6.2) 


The first derivative of the function is 


y (x) = ^ (-5y„ + 4yi ■ ^2^ * ? >^2) 

(6.3) 


The slope at each of the points x^ is obtained by substituting 
X = x^ = 0 , X = x^ = h, X = = 2 h in equation (6.3). The second 

derivative of the curve y(x) is 
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y"(x) = ^ (y^ - 2y| + y^) (6. 4) 

h ^ 

which is constant because y(x) is a second degree curve. 

Polynomials passing through foiir points and five points can be 
obtained in a similar fashion and their derivatives evaluated at each 
point to obtain various difference patterns. Thus, numerous choices 
are available when selecting a difference pattern. T'Jhich pattern is 
best depends to a large extent on the equation to be solved and its 
boundary conditions. However, one set of central difference operators 
is usually suggested in textbooks (refs. 4, 10, and l6) , widely used 
in the literature (refs. 2, 12, and l4) , and generally accepted as 
preferred because of simplicity, ease with which boundary conditions 
are handled, and consistency of order of error. These are given in 
equations (6.5) to (6.8) 

‘ S (-^1-1 

^ (^1-1 - * i^i+i) 

^ (*i) “ ^ (‘^1-2 * ^1-1 ° ' '^1+1 ''' ^i+a) 


(6.5) 


( 6 . 6 ) 


(6.7) 


y"''(q) = ^ (^i -2 - '*=' 1-1 *■ ^=*1 - ‘*=' 1+1 ■* =*1+2) 


( 6 . 8 ) 
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The simplicity of the selected operators equations (6.5) to (6.8) 
compared with writing at each station equations (6.3), (6.^), or those 
obtained from polynomials through four or five points (see ref. 16), 
to obtain the various difference operators is obvious. The order of 
the truncation error for each operator may be obtained by expanding 
the function yf about the point x^ in Taylor series as given 
below 
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where y^ stands for the derivative a is any real number and 

dx 

h is the increment of the interval. For example, equation (6.5) 
which is (6.5) evaluated at the center, can be expanded as follows 
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Che truncation error is of order h and is 


iii 

T^i 


V 

120 ^i ^ 


( 6 . 10 ) 


For the first difference, equation (6.3) evaluated at the left end 
point leads to 

4 (-5^1 *'^1+1 + ^1+2) 


and the truncation error is 


i£ iii 
5 ^1 
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(6 . 11 ) 


The second difference, equation ( 6 . 4 ) evaluated at the center point, 
is equation (6.6). It yields a truncation error of 


h^ iv 2l vi 
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( 6 . 12 ) 


The second difference, equation ( 6 . 4 ) evaluated at the left end point 
gives the truncation error 


hy: 
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(6.13) 


Note from equations (6.IO) to (6.13) that while the error for 


both first difference operators is of order ii , the error for the 


second difference operator about an end point is of order h, and 
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2 

about the center point, of order h . Generally, difference patterns 
with the same order of error are used for more consistent results. 

p 

For example, one of the first difference patterns (error of order h ), 
is not mixed with the second difference pattern at the end point 
(error of order h) . When inconsistent order of error terms are used, 
the answers tend toward the more inaccurate terms (ref. l6) . 

Difference Operations 

A second method for obtaining the various finite difference 
operators is differencing differences. The inverted delta, V, 
designates backward differences, the normal delta, A, forward 
differences and the lower case delta, 5, central differences. Suppose 
the values f^ = ^(^i) ^ function f(x) are known at (u + l) 

equidistant points = a + ih where i = 0, 1, 2, n 

(sometimes i is nonintegral). On the Interval (a, b) h is the 
Id 

increment and is taken to be positive. For any function f(x) 

the difference operators A, V, 6 are defined for increment h as 
follows 


^i = ^itl - "i 


(6.14) 


Vf. = f. - f_ 
1 1 i-1 


(6.15) 


" 1 

1 +^ 1- 2 


(6.16) 
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The differences in equations (6.l4) to (6.15) may be extended to 
higher order differences by taking the difference of the difference; 
for exang)le 



A /'-P 



2^i+l 


± . 
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(6.17) 


In general 


f^ = f^^ = v(vP"^ fj 6 ^ = 5(5^"^ f^) 

V ^ 1, 2 ■ ■ • 


(6.18) 


for p = 0 


A° f, = V° f. = 6° f. = f. (6.19) 

i 1 1 1 

Given in equations (6.20a) to (6.22e) are the finite difference 
operators that approximate the various order derivatives (including 
zero order) . Also included are the operators expanded in Taylor series 


to obtain truncation error terms . 
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Forward and backward differences give xmllateral expressions for 

the derivatives of a function > which in their simplest form 

have errors of order h. Central differences, involving pivotal 

points or stations symmetrically located with respect to x^, are 

parti c\jlarly useful in the solution of boundary value problems, 

reference 3- Note that the order of error for the central difference 
2 

operators is h . Generally, as h approaches zero the central 
differences approach the exact value faster than forward or backward 
differences. 

The central difference operators (6.22a) to (6.22e) are defined 
at half stations for odd derivatives . These operators are regular 
and consistent and may be used successfully in boundary value 
problems. They will be referred to as "half station" operators. 

The linear second order differential equation in the form 

L(y) = a^ (x)y" + a^ (x)y' + a^ (x)y = b(x) (6.25) 


cannot be approximated by the half station operators because the 


approximation for second derivative introduces unknowns y 


'i+1 


'i-1, ^1’ 

and the approximation for first derivative introduces unknowns 


y 2 ^ } y q • This gives too many unknowns for the number of 

i-— in — 

2 2 

equations. However, equation (6.25) nay tie reduced to the form 


L(y) = [^T(x)y'^ + g(x)y = p(x) 


(6.24) 
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Equation can be solved by using half* station operators. 

Averaging Procedure 

The central difference operators (6.5) to (6.8) which were 
obtained from Lagrange's interpolation formula, and which do not have 
half stations in the approximations for odd order derivatives, can be 
obtained by averaging difference operators. The first averaged or 
mean difference at i is obtained by taking the average of the first 
central difference at i--^ and i+^ . The operation is symbolized 
by the operator n called the averager. The first averaged 
difference is 


yj = I 
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(6.25) 
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Similarly the next averaged difference is 
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Since the averaged central difference operators are defined only at 
integral points, that is at whole stations, they shall be referred 
to as ’’whole station" operators. 

Expanding tht whole station operators in a Taylor series to 
obtain the first two truncation error terms results in the 
following 


Derivative Finite difference pattern Taylor series expansion 
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This study is concerned with determining and compaxing the 

accuracy of the two central difference methods both with order of 
2 

error h , the half and whole station methods. Also a modified form 
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of the whole station method, in which all derivatives occurring in 
the given equation L(y) are approximated, is considered. 
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VTI . A Method for Determining the Accuracy of the 
Central Finite Difference Equations 

To study the accuracy of the half station and whole station 
methods in second order boundary value problems, the simple problem 

L(y) = - gy = p(x) (7.1) 

with boundary conditions 

y(a) - y^ y(b) = y^ 

on the interval fa,b) is considered. 

13“ £l 

With h = " and = a + ih the finite difference method, 

applying operator equation (6.l6) to equation (T.l), and noting that 

-M'i 'K- M, 1+ M. r 

- 7 1 5 (- i'l-i 7 1 j (- 3^1 + yi+i) 

yields 

1. Half Station Method 

+ Sl^l - Pi - 0 

(7.2) 

2, • • • n - l) 
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Using the expanded form of equation (7-l), 

L(y) = - f'y" - fy' + gy = p(x) (7.5) 


and the operators in equations (6.5), (6.6), the difference equation 
for (7»l) takes the form 

2a. Whole Station Method 



2^i^i 




Si^i 



= y„ 



(7A) 

(i = 1 , 2 , • • • n - 1 ) 


The derivative f^ in equation (7*^) can be evaluated exactly at 
the appropriate stations. Another method which can be considered is 

I 

to approximate f^ by equation (6.26b). The result is 
2b . Modified Whole Station Method 



+ %yi - Pi = 0 (7.5) 
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Note that the three sets of finite difference equations, (7.2), 
(7.4), and (7*5)» lead to different coefficients for the simultaneous 
equations in terms of the same displacements at the i^^ station. 

With a few assumptions the existence and uniqueness of the solution of 
each of the sets of simultaneous eqAiations is established from a 
theorem proved by Collatz and stated in appendix A. 

If it is assumed that f(x) >0 and g(x) ^ 0, the systems of 
equations (7.2) and (7-5) satisfy the conditions of Theorem A1 (in 
addition to the sign distribution, the weak row sum criterion is 
satisfied and the matrix of coefficients is irreducible); hence, a 
uniquely determined solution exists for each system for arbitrary 
boundary conditions and arbitrary values of p^^. For the set of 
equations (7.4) the additional assAunption that for 


f'(x) 4 0, 



satisfies the conditions. 

The usual approach in a finite difference accuracy study 
(ref. 12) is to carry out the numerical solution to a number of 
problems for which the exact solutions can be obtained and compare 
the resulting numerical answers with the exact answers. This 
procedure was carried out for a number of problems of the type of 
equation (7.1) and a table of relative error for a typical resilt is 
given in appendix B. Such a procedure has the liability that cal- 
culations must be redone each time the increment size, h, changes. 
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The root mean square of relative error for three different values of 

the increment h for problems solved is given in appendix B. 

Conventional means for estimating the error bounds do not give a 

satisfactory means of comparison of the different methods since the 

error limits exceed the actual error in magnitude. 

To obtain definitive expressions for error in each method, 

independent of increment h, first expand the finite difference 

recursion equations (7.2), and (7.5) in a Taylor series 

"til 

expansion about the i point. For each method this leads to a 
differential equation of the form 


^o (^i) - Pi + ^ ^ ^2(^1) + • • — 0 (T.6) 

subject to the boundary conditions 



^i = ^a 

at 

X = a 



^i = ^b 

at 

X = b 


The symbols 

^0^ i"!# ^2 

are linear 

differential 

operators given 


^o(h) 

= - (vi' 

)' % h 

(7.7) 


and 


1. Half Station Method 
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I-V -*-/ 


iv iii ,, 

Vi . Vi ^1 . ^i ^i 


(7. 8a) 


^2(^i) 


I l“i ^ i^i ^ 'i^i ^ i ^i ^ i ^i ^ li^i ) 

1 5^0 120 96 lU 58^ 1920 j 


2a. Whole Station Method 




iv iii 

__ + __g — 


\ 


(T.8b) 




" I 360 


fJy^ 

i’'i 

120 


2h. Modified Whole Station Method 




11 . 11 


N 


,12 6 

^ J 


(7.80) 




- vi X., V „iii iii _v , 

Vi , ^l^i , ^i ^i , VI 

3& 120 120 56 


Equation (7*6) and (7.7) together with (7.8a), (7.8b), or (7«8c) 
are clearly the differential equation which represent exactly the 
finite difference equations. As h approaches zero, equation (7-8) 
approaches equation (7.I). The solution to equation (7.6), satisfying 


- 50 - 


the appropriate boundary conditions, gives an analytical representation 

of the numerical finite difference answers. A closed form solution to 

equation (7-6) does not appear feasible since it contains an infinite 

number of terms. For a practical problem, however, if the length of 

the interval (a,b) is one unit, h is perhaps 0.1 or 0-01 or even 

smaller. This suggests that equation (7-6) can be solved with the 

2 

use of perturbations with the parameter taken to be h . 

Let the solution y^ to equation (7-6) be taken in the form 

y^ = + h^^ + . . . (7-9) 

Substituting equation (7-9) into equation (7-6) leads to 



Pi + 






+ 


• = 0 (7.10) 


subject to 


Y^(a) + 

Y^(a) 
_ — 

Yo(b) + h^ 1 

Y^(b)"^ 


If each order of error term is solved in sequence, the following 
series of problems result - 

Y^(a) 0, Y^(b) - 0 




Pi = ° 


(7.11) 
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(2) L^(Yi) + L^(Yj = 0 Y^(a) = 0, Y^(b) = 0 (7.12) 

( 5 ) • • • 


Note that since eqiiation (T-l) is linear Y^ given by 
equation ( 7 .II) is in fact the exact solution. From the form of 

it is seen that Y^ can be interpreted as the first order error 
term in the finite difference results. The magnitude of Y^ is, 
therefore, a measure of the error in the finite difference results as 
compared to the exact answer to the problem. A comparison of the 
error terms Y^ resiilting from the different finite difference 
approximations indicates the relative accuracy of the different 
approximations . 

While errors in the y^ are important, errors in numerically 
obtained derivatives should also be considered for a thorough error 
analysis. Therefore, results were obtained by using the finite 
difference answers for approximate second derivatives. The second 
difference operator was applied to the difference results followed 
by Taylor and pertttrbation series expansions to yield 
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VTII . Application of the Method to Particular Problems 

Problems Studied 

Using the method described in the previous section, the error 

term in equation (7*9) and the second derivative error term 

yiv 

Y'-^ + in equation (7.15) have been obtained for a series of 

problems for the half station and whole station approximations. 
Equation (7-l) has been solved with g = 0, p - - 1 for the 
following values of f(x) 


(1) f(x) = ^ for 

X 


subject to the boundary conditions 


l<n:^6 i5x<2 

( 8 . 1 ) 


y(l) - 0 

y(2) - 0 


and 


(2) f(x) = 1 + x“ for 


subject to the boundary conditions 


l^n£5 O^x— 1 


(8.2) 


y(0) = 0 


y(i) - 0 
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Physically these problems might correspond to the problem of lateral 
deflection of a string having a uniformly distributed lateral load 
and a variable tension force f(x) . 

For the case where f(x) is linear (corresponding to f(x) = 1, 

X, or 1 + x) the results for the half station and two whole station 

finite difference approximations are exactly the same. In fact for 

f(x) =1 1, all three difference answers are the exact answer. For all 

other cases, however, the three difference methods lead to different 

results. It is useful to compare the results for the case 

f(x) = in detail as a typical example, 

x^ 

For f(x) = ~ and y(l) = y(2) = 0 


Y 

o 


5 75 “75 


( 8 . 5 ) 


and 


1 . Half Station Method 


^1 4 x^ 11 2 , 86 

1 “ ■ 1125 * "S" " 150 * 1125 


( 8 . 4 ) 


2a. Whole Station Method 


^1 " 5 50 225 


2b. Modified Whole Station Method 


Y _ _ 101 x^ + 

1 “ 90 5 ^ 55 


( 8 . 5 ) 
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A plot of the three error terms over the vmit interval is 

given in figure 2(b). The exact solution, is given in figure 2(a), 
The finite difference solution (7-9) can be obtained to the first two 
terms for any desired increment h from figures 2(a) and 2(b) . 
Solutions were also obtained for the error terms for i of the 
remaining fimctions f(x) noted previously) additional plots of 
results and the exact solution for the case f(x) = 1 + x^, are 
shown in figures 3(a) and 5(b). Detailed plots of the remaining 
solutions are not shown because figures 2(b) and 3(b) serve to 
illustrate the character of the resiiltsj and overall measixre of the 
relative errors in the two methods will be shown for all the 
solutions obtained. 

The error terms for the second derivatives corresponding to the 
different methods and for the case f(x) = are as follows 

1. Half Station Method 


Y” 



l64 2 
575 ^ 


X + 


51 

75 


( 8 . 6 ) 


2a. Whole Station Method 


Yiv 

Y" + — 
1 12 


^ x2 b 6x - 


75 


25 


(8.7) 


2b. 


Modified Whole Station Method 

,lv 


Y 

Y^ + — 


206 2 ^ 

— X + 26x - 

15 


715 

75 


12 


( 8 . 8 ) 
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A plot of the error in the second derivative for each of the methods 
is given in figure 2(d) for this case and in figure 3(d) for the case 
f(x) = 1 + x^. The exact solutions are given in figures 2(c) 

and 5(c)- Again the first tvo terms of the finite difference 
solution (7-13) can be obtained from these plots for the desired 
increment h. Results for the remaining functions will be shown 
later . 

Numerical calculations were also carried out for the deflections 
and the second derivatives for the problems cited to determine if the 
analytical errors adequately represented the numerical errors. The 
data are not included herej however, for h less than about 0.1 all 
numerical errors agree with analytical errors to within one percent. 

Relative Errors of the Half and Whole Stations Methods 

While resxilts such as those given in figures 2 and 3 are usually 
sufficient to identify which of the methods is superior for a given 
problem, identification of the superior method for specific results is 
sometimes difficult. A quantitative measure of the relative accuracy 
of the methods can be made by examining the root mean square values of 
the errors for the entire solution, that is 



for the error in deflection and 
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for the error in second derivative, where the integration is over the 
unit length from a to b. Thus, to assess quantitatively the relative 
merits of the half station and whole station methods for the various 
problems solved, the ratios 

^l,half 

^1, whole 

Y*' 

l,half 

T' 

1, whole 

have been calculated for each problem. The results are shown in 
figure 4, Ratios for the modified form of the whole station method 
compared with the half station method are given in figvire 5, 

Discussion of Results 

The results given in figures 4 and 5 show that for the problems 
studied the error in the deflection resulting from use of the half 
station method is less than the error resulting from the use of the 
whole station method. The investigation of the accuracy of the 
second derivative approximations gives the same result in general . 
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The difference between the two methods is generally less in 
calculating the second derivatives of deflections than in calculating 
the deflections themselves; moreover, differences in the comparative 
error from problem to problem are noticeably less with the second 
d,62!*xva.‘tiy6S 'tlifi.n wx'tii *t}i6 d©n.0C'bioiis* 

It should be noted that the analytical representation of errors 
shows clearly the danger of using numerical data at a single station 
or a few points to characterize the error in a problem. A typical 
case is shown in figure 2(d) for f(x) = . If comparisons are made 

of the second derivatives near the end x = 1, the whole station 
method appears much more accurate than the half station method; 
however, figure 4(b) shows clearly that the average error with the 
whole station method is more than twice as great. 

Reasons for the superiority of the half station method are not 
altogether clear, but may Include the symmetry of the matrix of 
coefficients in this method. By contrast, the matrix of coefficients 
associated with whole stations is not symmetric. Matrix symmetry can 
be of great value for many nmerical procedures associated with 
eigenvalue routines and simultaneous equation solving routines and, 
in some cases, is required for an efficient numerical solution of a 


large order system. 
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IX. CONCLUSION 

A procedure was developed to determine an analytical expression 
for the discretization error in a finite difference solution to allow 
a direct comparison of methods which was independent of the increment, 
or mesh size. Using this procediire, a comparison was made of the 
accuracy of two different finite difference methods for solving 
linear second order boundary value problems . 

The methods investigated were a "half station" method which 
corresponds to making the finite difference approximation before 
expanding the derivatives of function products and a "whole station" 
method which corresponds to expanding such products before making the 
approximations. Both of these methods are currently in use. Also 
investigated was an alternate form of the whole station method in 
which known derivatives are approximated rather than evaluated 
exactly. It was found that, for the same number of stations, the 
average error in calculated deflection resulting from use of half 
station difference approximtions was always less than the error 
which resulted from the use of the whole station difference 
approximations. In some cases this error is reduced by an order of 
magnitude. It was also found that the alternate form of the whole 
station method gave the same or better results than the usual whole 
station approximation. The investigation of the accuracy of second 
derivatives gave similar results in general. 
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X. BIBLIOGRAPHY OF SOME PUBLICATIONS ON NUMERICAL METHODS 

Numerous publications concerned with the numerical solution of 
differential equations and the accuracy of these solutions are foxind 
in the literature. Those publications that are useful in the present 
study of the accuracy of central finite difference methods for 
approximating the solution of boundary value problems in ordinary 
differential equations are listed as references in this paper. In 
addition a number of publications which are concerned with the numerical 
solution of initial value problems, or problems which can be changed to 
this type, and the numerical solution of partial differential equations 
are given in the bibliography. 

The bibliography is arranged in five sections. Included in the 
first section are publications in which the theory of one or more of 
the different methods for obtaining numerical solutions is discussed. 

In some of these articles, discussions on stability, convergence and 
accuracy are included. The second section includes publications in 
which the emphasis is placed on error estimates and error bounds. The 
third section contains publications which report on the methods used 
to obtain approximate solutions of ijarticular physical problems. In 
the fourth section are books on methods of numerical analysis. The 
last section contains some extensive bibliographies which cover the 
different areas of numerical analysis. 
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XIV . APPENDIX A 


The following theorem is proved by Collatz (ref. 5? page 44) . 
Theorem A 1. If the coefficients a , of an n x n matrix A 
satisfy the conditions 

1. Sign distribution a >0, a < 0 for j / k 

J J 

2a. The weak row- sum criterion 


^ a 

L 

^ 0 

for 

.3 = 1, 2 . 

k=l 

> 0 

for 

at least one 


and 2b. Matrix A is irreducible or instead of 2a and b the 
stronger condition 

2c. Ordinary row sum criterion 


n 



k=l 


then A is monotonic and det A ^ 0. Thus a unique solution to 


A exists . 
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XV. APPENDIX B 

The problems named in equations (8.1) and (8.2) were solved 
exactly and numerically by the half station, whole station, and 
modified whole station methods with the increment h eqxial to 0 . 25 , 
0.125, and 0 . 0625 . The relative error R was calculated for each 
method- An example of the table of exact and approximate solutions 
and the relative errors at several points is given in table I for 
the problem -((l + x^)y') = 1. Rather than include similar tables 

for each problem, the root mean square of the relative error given by 

■ ■ ■ *4 

was found for each method. The root mean square errors for each 
method for each increment h are given in table II. After examining 
tables I and II it is seen that it is sometimes difficult to 
determine which method is best for the desired increment h. 


EXACT SOLUTION, APPROXIMATE SOLUTIONS AND RELATIVE ERRORS OF APPROXIMATE SOLUTIONS TO THE 


* 


6l 


o 

II 


>> 

II 



r-i 

II 


>> 


X 


+ 

‘H 


I 


B 

H 

I 


•0 




fi 

<C <0 O' 

'C 0* ^ C <M 


h V 

® r* •- 

<C >C <0 <C ^ IT 

^ ^ «p4 ^ 


• • • 

• ••»••• 


C O 5 

fVj <V CN. 

o o c coco 

6 c c o c o o 6 c o c o o C' o 

w a ? 
>(»- 




u « 

tsi 

O OD iT» OC' 00 0^ 


So 

•t <D 

tr IT* ^ <c <o IT. 

^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ «m4 


• • • 

• •«•••• 



rg <NJ fsj 

C O O C O CO 

coobcocodoooooc 





u 

O' OC IT 

C ^ ^ ^ 

irif. 'C'O'CCP-O'C'C^U'iT'd-'i' 

8 H 

a C 00 

r»y (V! rs» r*j rg CN* ^ 

cccccooccccocco 


• • • 

• ••••«• 

• •••••••••••••• 

H Js 

c c 

o c o o o c o 

coocooccocooooo 





s a 




55 

•& O' ir\ 

^ c ^ ^ ^ 

O— ‘O'tir— 'Cif'— i^C'j-cco'if'. 

•rJ 

'S‘§ 

^ ec 

O If* O' IT* ^ ® 00 

ocBcspccO"-‘ooir«c--rj'£)-4->4-PJ 

O' 

1^ «c i*> •“ fvi ^ n 


ir or o 

^ ^ yr* ^ If 

fv cv IT. c IT r~ >c 00 ec O' •«■ ■♦ OC 

0> Q " 

IT 

cc ^ ^ ^ 

c «- Pw ^ X' p- ^ o c p- IT c If c 

•H X 

fsj C 

CD C 00 ^ 

^flC'£>CC^0'P'r''XXIO'0'fPP* 

Vi w o 

•iH t— ! M 

or o »c 

^ tc O- O' ct *c 

(V •» «C' OC o O' O' O' 0 OC p- X fp t-< 

c — < c 

o c c o c c c 

000000000000000 

0 5 S 


• •••••« 

• •••••••••••••# 

g? 5* 

o o c 

o o o c o c c 

ccooocccocoooco 

a a 
o o 

^ < 

«C f»'i ^ 

o^rr. o Xl<v.'i'fMO»-> 0 ''C> 0 ' 0 'C 


<Si ir •- 

IT IT O tf' >t a 



ic < 

^ o ^ >c >c 

PvJfVX— 'C'i’O'CO'OC— 

43 § 

O «'■ O' 

C fT K- <4- fr ^ C 

.'-CP’‘Oc^^xxf'-0'0--xxO 

m <H 

fVJ U" 

<NI <Si ^ <Ni IT .^- 

o PO <*■ X f>- O' rr, p- X fv O X O 

X 

fg O 't 

cr o ec ^ f»^ 

*C«%CCO'CO'f^frx. Iff^o 

3 8 

or. o < 

•tf- X* 0“ O' X 

fV ^ »C X 0 0 O' O' 0 X %C ij- ff ^ 

O ft 

c ^ c 

o o c o c c c . 

COOCOCOCCCOCOCO 

5 ft 


• •••••• 

• •••••••••••••• 

5 ff 

c c o 

o c o e c c c 

COC'CCCCCOCOOOCO 

rt S 


IT X. IT ff ^ X ^ 

X fv O CV OC X' fvi J- X O' X. O »- (V fp. 

^ P 

fT 

^ C X ^ ® X *c 

x»-pvipv — p..crPO xc^pp. X. •-< 


or u" 

f i X ^ (V fw ^ 

X' X X X a (V c o c X c X x pv ^ 

4i oS 


^ tr ic ^ <c «*• >c 

O' «c X — X. p- f' "S' X o fv o fv 1" 

os a 

O CC ' 'C 

-< S' c op ^ ir 

15' c ^ fT. ^ ^ C X fv '^- C X o 

•*3 -H 

^ tt ff. 

a C ff Xt rr 

X. X X' c c >o o p“ X X. X x> o X 

CD X 

X a ^ 

.S' X 0 0^ X ^ rf 

pvX-xxc'O'O'O'cxp-X'tx. •- 

Vi S 

c c c 

c c c c c c 

C C C C C C C C C O C C C c c 

r-l 


• •••••• 


ss* 

o o o 

o o o o o o o 

OOOOOOOOOOOOOOO 


f<v 

O' rCL f^^ if\ c«% O 

X0''4'XiXXPJXXX<0X— 400 


4- .JS fw 

-0 ^ an O' (V tr 

ob^'#'xcoo'0'-4fvi>j'p'XX30 

s 

iT^ tr» o 

If 0- if if O X 

■^XOXPOO'XXPsIXI^OXCOf'- 


g; oD -« 

-t -O O so IT' -- o 

c04'-4'OOOO®PslXP-»4a0OO 


fN» ^ 

O tM >0 •■* If 

0'0'^<M'l'OOP»P'i0 4'»^0'XO 

O N. 

» O >0 if ^ 

xaooooo'f*-xxxx®xp- 

rt 

TD ^ «0 

X ^ 0. CO ^ ff 

P4-tO»0'0'^3'O'®P“O'*'X-^ 

^ o 

c c c 

C C O C C C C: 

CCOCCOCCCC'OOCCO 


• • • 

• •••••• 

• •••••••#•••••• 


c c c 

C C C‘ C c o c 

c c c c c c c c c c c c o c o 


r. c c 

c c c c o c c 

xcxcxcxcxcxcxcx 


c c c 

X r IT c X r IT 

X X c X' X p- c XI X P- c X‘ X 


IT. C U* 

PJ IT, p- c C-' ir 

xxxx '-'t' xe>oxxx*-'X-x 

X 

rv ir 

•-» f ff“ X 'f X 

C ^ X X X -f X X O O X X OC O' 


o o c 

c C C. c c o c 

c c c c 6 b c o o o o c o c 



- 62 - 


TABLE II 


ROOT MEAN SQUARES OF THE RELATIVE ERRORS OBTAINED 
IN APEROXIMATING THE EQUATION - (f(x)y')’ = 1 


f(x) = i 

X 


Size of increment, h 

1/4 

/o 

I/O 

1/16 

% error using half station method 

.42 

.16 

.06 

^ error using whole station method 


.65 

.24 

^ error using modified whole station method 

1R1 

.46 

.16 


f(x) = ^ 

X 


Size of increment, h 

1/4 

1/8 

1/16 

^ error using half station method 

.51 

.12 

.04 

^ error using whole station method 

5.82 

1.43 

.52 

% error using modified whole station method 

.74 

.74 

.28 


f(x) = 


Size of increment, h 

1/4 

1/8 

1/16 

^ error using half station method 

.59 

.16 

.06 

% error using whole station method 

5.91 

2.15 

.77 

^ error using modified whole station method 

5.95 

1-95 

.78 


f(x) = ^ 

X 


Size of increment, h 

1/4 

1/8 

1/16 

% error using half station method 

1.62 

.64 

.24 

% error using whole station method 

7.45 

2.57 

.92 

% error using modified whole station method 

12.00 

5.52 

2.17 































































- 6 


T 


f(x) ^ 


r- — 

Size of increment, h 

lA 

1/8 

1/16 

% error using half station method 

3.^0 

1.54 

.50 

% error using whole station method 

8.15 

2.79 

1.00 

/S error using modified whole station method 

28.22 

12.10 

4.69 


f(x) = 

^ ' D 

X 


Size of increment, h 

1/4 

1/8 

1/16 

^ error using half station method 


2.24 

.85 

io error using whole station method 


5.58 

1.28 

io error using modified whole station method 

52.40 
1 

22.11 

8.55 


f(x) 1 + X 


Size of increment, h 

1/4 

1/8 

BE3' 

% error using half station method 

.90 

.35 

BBII 

% error using whole station method 

.90 

.55 

.14 

% error using modified whole station method 

.90 

.55 

.14 


f(x) = 1 + X^ 


Size of increment, h 

1/4 

1/8 

1/16 

% error using half station method 

1.16 

.45 

.17 

% error using whole station method 

3.25 

1.25 

.45 

% error using modified whole station method 

5.25 

1.25 

.45 


f(x) = 1 + 


p— — — 1 

Size of increment, h 

1/4 

1/8 

1/16 

io error using half station method 

1.69 

.5^ 

.22 

ia error using whole station method 

4.58 

1.64 

.58 

% error using modified whole station method 

4.54 

1.62 

.59 


TABLE II.- Continued 
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k 

f(x) 1 + X 


Size of increment, h 

lA 

1/8 

l/l6 

^ error using half station method 

2.18 

.78 

.28 

error using whole station method 

5.79 

2.08 

•75 

error using modified whole station method 

k.ik 

1.35 

.47 


f(x) = 1 + x^ 


Size of increment, h 

1/4 

1/8 

1/16 

% error using half station method 

2.72 

.97 

.34 

^ error using whole station method 

7.12 

2.56 

.92 

^ error using modified whole station method 

3.05 

.88 

.31 


TABLE II.- Concluded 
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Figure 1.- Function y = f(x) . 
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Figure 4 .- Ratio of root mean square errors for half and whole station methods. 
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Figure 5»- Ratio of root mean square errors for half and modified whole 

station methods. 



